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Abstract. We construct and analyze a preconditioner of the linear elastiity system dis- 
cretized by conforming linear finite elements in the framework of the auxiliary space method. 
The auxiliary space preconditioner is based on discretization of a scalar elliptic equation with 
Generalized Finite Element Method (GFEM). 



1. Introduction 

The discretizations and fast solvers for the linear elasticity systems have been extensive 
subject of research for the past 20 years. A number of methods have been introduced 
in [3l El [121 E31 EH E], and in these works one can find stable discretizations, a priori error 
estimates, as well as construction of fast and robust solvers for different range of material 
parameters. 

In this paper we consider the variational problem corresponding to the lowest order finite 
element method discretization of linear elasticity system in displacement formulation. We 
assume that the material parameters in the linear elasiticity system are well behaved, namely, 
the Poisson ratio is away from 1/2. We will not discuss the robustness of our method with 
respect to these parameters, because our goal here is to introduce and prove results on 
the relationship between GFEM discretizations of scalar equations and the linear elastitcity 
system, as well as to employ such relations in the construction of auxiliary space [151 f!8] 
preconditioner. 

The preconditioner that we construct employs as auxiliary space the piecewise quadratic 
conforming finite elements, and the corresponding auxiliary bilinear form on the auxiliary 
space corresponds to the discretized scalar Laplace's equation. To relate the discretization 
of the scalar Laplace problem to linear elasticity, we also use an "intermediate" GFEM space 
containing piece- wise quadratic functions. The key steps in the analysis of the preconditioner 
rely on the spectral equivalence results proved in [11] and the fact that the kernel of the 
quadratic GFEM stiffness matrix is isomorphic to the space of rigid modes (see [TT]). 

The remainder of the paper is organized as follows. In section 2, we present the linear 
elasticity problem of interest and its variational formulation. In section 3, we show an 
auxiliary spectral equivalence relation that plays a key role in the analysis. We introduce 
the preconditioner, and prove uniform spectral bounds in section 4. Numerical results are 
presented in section 5 are shown to validate the theoretical results. 
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2. Linear Elasticity 

Let Q C IR 2 be a bounded convex polygonal domain with boundary V = dQ. We consider 
the linear elasticity problem with pure traction boundary conditions: 

-div(2/ie(u) + Xtr(e(u))S) = f in fi, 

(1) . . 

(2fie(u) + \tr(e(u))5)n = g on T, 

where / is an external force, n is the outward unit normal on the boundary, = ^{djU{ + 

diUj) is the strain tensor, and 5 is a matrix whose elements consist of the Kronecker delta 

symbol. Further, let tr denote the trace of a matrix and define the Lame coefficients /i and 
A in terms of the Young modulus E and the Poisson ratio v as follows: 

A = — and /i 



(1 + i/)(1-2i/) ^ 2(l + i/) 

We mention that due to the corners of the boundary of the polygonal domain Q, care must 
be taken when considering the boundary conditions of ([I]) (See [8] for details). 

Let Si, 1 < i < n, be the vertices of r, Tj, 1 < i < n,the open line segments joining Si to 
Si+i, and ni the unit outer normal along Tj. Let p E H 1 ^ 2 ^^ and q E H l l 2 {Yi + i). Then, 

p = q at iSi + i if 



f 5 ds 

/ |?(«) — P(— s)| 2 — < 00, 



where s is the oriented arc length measured from Si+i, and 5 is a positive number less than 
min{|rj| : 1 < i < n}. Then, equation (CQ) can be written more precisely as 

-div(2/ie(w) + Atr(e(u))|) = / in 

(2) 

(2/ie(u) + Atr(e(tt))5j7}i| r = ^ 1 < z < n, 

where / E L 2 {Q) and gi E H l/2 (Ti) satisfy 

g^i • n i+ i = g i+ i ■ Hi at S i+1 for 1 < % < n. 

Denote by RM and RM -1 the space of rigid modes and its orthogonal complement with 
respect to L 2 (f2), respectively and by H]_(£l) := H 1 ^) flRM -1 the intersection of H^Q) with 

the L 2 orthogonal complement HM ± of RM. From integration by parts, we obtain the weak 

formulation of (J2J) as follows: 
Find u E H]_(£l) such that 



/ / ■ v dx + / gi- v\ Tl ds 



(3) a LE(u,v) 
for all v E H]_(Q), where 

a LE (u,v) := J {2^{e(u),e(v)) F + A(V • u)(V ■ vfj dx. 
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Here and below, (-,-)f and ■) are the Frobenius inner product for matrices and the standard 
Euclidean inner product for vectors in M. d , respectively. The corresponding norms are denoted 
with | • \p and | • |, respectively. A sufficient condition for existence of a solution to ((3]) is 
that the following compatibility condition is satisfied: 

J f-vdx + 'Y^J g t ■ v\ r ds = 0, Vt> G RM. 

Following [19], we write X\ < y\,x 2 > y 2 and x 3 ^5 y 3 whenever there exist constants 
Ci, c 2 , c 3 and C4, independent of the mesh size h, such that 

xi < Cij/i, x 2 > c 2 y 2 , and c 3 x 3 < y 3 < C 4 x 3 . 

Moreover, we write <a? and ^a when <, > and ^ are dependent on the the Lame coeffi- 
cient A. 



3. A MULTILEVEL PRECONDITIONING FOR THE GFEM PROBLEMS 

In this section, we present our GFEM-based preconditioner, using results from [H] to 
motivate our scheme. Let Q C lR d be a polygonal domain, with d = 2, 3 and V\ be the 
associated piecewise linear finite element space on a quasi-uniform triangulation Th of Q. 
Then the GFEM space V GFEM is defined as 

yGFEM := Vl + J2Yl WlMi}> 
k=l i=l 

where V\ corresponds to continuous piecewise linear elements {0j}" =1 and ip k ,i — Xk h k ' 1 with 
Xk the fc-th component of x and x k) i the k-th value associated with the i-th nodal point. 
Consider the bilinear form 

(4) a(u, v)= [ Vu ■ Vv, Vu, v G V GFEM . 

Jn 

It is shown in [TT] that the kernel of the stiffness matrix A corresponding to (pEj) is charac- 
terized by the space RM of rigid modes. To obtain an efficient preconditioner for a(-, ■) on 
yGFEM ^ we j-^g ^QQgg as auxiliary space W = V\ x Vf defined by 



n n 



u G W, u =(u,u), u=(u k ) where u = }] a^j, u k = y^ j a k ,i<j>h 

i=i i=i 

equipped with a w (u, v) = a(u,v) + J Q (e(u), e(y)) F dx. Define II : W — >■ V GFEM to be: 

n d n 

u -> flu = ^ + ^ a k ,iilJk,i<fii, 



t=l fe=l i=l 



and take A G : \/ GF£:M _>. (yGF^My anci ^ : W' to be the isomorphisms associated 

with a(-,-) and a^(-, ■), respectively. Here, ' and * denote the dual spaces and adjoint 
operators. 
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We consider the case when the linear part of the GFEM elements are zero: W = {0} x Vf, 
so that II : W — > V , or more precisely, 



u = (o, u) — > riu = o + ^2 ^2 ak ^> 



d n 



k=l i=l 

or 

d n 

uew ^Uu = J2J2 a ^Mi e v GFEM . 

k=l i=l 

We remark that ker(AG) is isomorphic to ker(A-w) where ker denotes the kernel of the 
operator. The following spectral equivalence, proved in the next subsection, then holds: 

(5) a^(u, u) ^3 a(ITu, Tin), \/u e W. 

Now, by the auxiliary space lemma [T5l HE], ni?n* is a uniform preconditioner for Aq, 
namely, k(UBU* Aq) < 1 where B is a BPX preconditioner for Ay/ (see [TT] for a detailed 
proof). 

A w 



w w 



II 



n* 



yGFEM y fyGFEM\l 

A G 

Hereafter, we use Hg for the GFEM interpolant associated with the GFEM space, V 



G 



3.1. Spectral Equivalence. Consider the following bilinear form defined for piecewise lin- 
ear, continuous vector fields on a triangulation of Q with simplexes Th 

a w (u,v) := / (e(u),e(v)) F . 



For a given piecewise linear continuous vector field, u, denote the corresponding element 

in the GFEM space by uq- Consider the isomorphism, II, between the GFEM space (with 
linear part) and the piece-wise linear continuous vector fields in W. d : 

d n 

u G = uu = ^2^2 a k,iipk,i4>i- 

k=l i=l 

By direct computation 

(x,n) - (x,n)j 

«g(x) = T , 

where wi denotes the continuous linear interpolant of a function w. We now arrive at the 
following lemma. 

Lemma 3.1. The following relations hold for any fixed simplex T of the triangulation and 
any x G T: 

(6) u G (x) = -((£(w)x,x) - (e(u)x,x)j). 
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Proof. The proof follows by taking a Taylor expansion of u and using the fact that u is a 
linear vector field on T and hence Vw is a constant matrix on T: 

(n(x),x) = (u(x ),x) + ([Vu]x,x) - ([Vu]x ,x). 

Taking the linear interpolant on both sides and subtracting leads to 

u G (x) = -«[V«]x,x) - ([V«]x,x) r ). 

The proof is concluded by observing that for any matrix Z e M dxd we have 

(Zx,x) = i((Z + Z r )x,x), 

where Z T is the transpose of Z, □ 

Writing out 1 = J2i=i ^« an d x = X)t=i ^i x i; where Xj are the vertices of T, and Aj(x) are 
the barycentric coordinate functions, we obtain that 

(7) u G = ^2(e(u) Xl Xj ,xt- Xj)tp E (x) 

EcT 

where <p E = K^j on the edge XjX]. Differentiating and taking the L 2 norm we have 

(8) ||V« ||g >T = J2 E { £ ^T> XE )([ V^"V^)(e(n)^,^ 

E'CTECT ^ T 

where The proof of the spectral equivalence lemma uses the following technical 

result. 

Proposition 3.2. Let Z G be a symmetric matrix and T be a nondegenerate simplex 

in M. d whose edges of size (independent on h). Then 

\ z \f ~ J2((ZyE,yE)) 2 , 

ECT 

where the constants of equivalence depend on the spatial dimension d, and yE is a vector 
parallel to the edge E. 

Proof. We first prove that the following expression is an inner product on Mfy^: 

(Y,Z).:=Y,{(ZVE,VE)(Yy E ,y E )). 

ECT 

First, we aim to establish that 

(Z,Z)*>0, and {Z,Z)* = iff Z = 0. 

We now provide a detailed proof for d = 3 case. A similar proof applies to the d = 2 case. 

Assume that (Z, Z)* = 0. This implies that (Zy E ,yE) = for each E C T. Since the size 
of T is independent of h, we can map it onto the canonical simplex T in WL d (the convex hull 
of the canonical coordinate vectors in IR d ). Clearly such mapping is affine and independent 
of h as well. Thus, we can limit our proof to T. Note that three of the edges of the canonical 
simplex are parallel to the coordinate vectors {ek}l =1 , and there are three more, which are 
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parallel to the vectors = (e^ + ej), 1 < k < I < 3. Let Z = I a\ 2 a 2 a 23 . From 

\a-i3 a 2 3 a 3 J 

(Zek, ek) = it follows that = 0. Next, {Zeki, &ki) = implies that a*,. + ai + 2(2^ = 0, 
and hence a&j = as well. Proofs of the remaining properties of the inner product are 
straightforward and hence omitted. 

The proof of the proposition is then concluded by using the fact that all norms on the 
finite dimensional space Mf*^ are equivalent, with equivalence constants depending on the 
dimension. □ 

We now prove that the Poisson bilinear form on GFEM and aw{-, •) are equivalent. 

Lemma 3.3. The following equivalence relation holds 

||Vu G ||o,n ~ ||e(«)||o,n- 

Proof. Take T G Th- Note that uq is a quadratic function on T, namely 

Ug = ^ ®e<Pe, where Q E = (e(u)^-, x E j , 

EcT 

This function vanishes at all vertices of T, and its linear interpolant (uq)i is zero. Thus, 

h ~ 2 \\ u G\\l iT < ll Vu G|lo,T ~ h ~ 2 \\ u G\\l, T i 

where the first inequality follows from the standard interpolation error estimate and the 
second is an inverse inequality. On the other hand 

II^gIIo.t ~ 1^1 E 

ECT 

Hence, 

(9) \Vu G \l * h~*\T\ £ 6| = |T| E (s(u)^, X f)\ 

EcT EcT 

Note that the size of ^ is independent of h and thus we can apply Proposition 13.21 (with T 
homotetic to T, and with edges y E = nf). Hence, from Proposition 13.21 and © we obtain 

IK«)llo,n= E IK«)ll2,r= E / <*(«)>*(«)>*■ 

(10) ^piE«t)'"Ewi 

j 

l|Vu G ||an. 



T£Tk EcT T&T h 



□ 

4. Preconditioning for Linear Elasticity 

In this section, we develop and analyze an efficient preconditioner for the linear elasticity 
problem. An additional spectral equivalence is needed to verify the auxiliary space lemma 
and, hence, the optimality for our preconditioner in this setting. This spectral equivalence 
result follows from a Korn inequality on H]_(Q) and the Cauchy-Schwarz inequality. 
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Proposition 4.1 (Korn inequality). There exists a positive constant C such that 

Mv)\\^>C\\v\\ 1>n We 

Proof. The proof is given in [H [7] . □ 

It follows from the Korn inequality that the weak formulation ([3]) has a unique solution 
u G H\(Q). 

Theorem 4.2. Further, the following spectral equivalence holds: 

ll*fe)llin S ^le(v,v) < a \\e(v)\\l >n , \/v G H\(Q). 

Proof. Let v G H]_(Q). Then we have 

a £i?fe V) = M\ £ (v)\\o,a + A|| V • u||g j0 

<MHv)\\l u + x\\v\\l a 

< 2n\\e(v)\\* fl +X\\e(v)\\l n 



<A Hv)\\l n , 

where the Cauchy-Schwarz inequality and the Korn inequality are used. The reverse inequal- 
ity is verified as follows: 

Mu)\\o,a £ MHv.)\\ 2 0>n + A||V ■ u||g )0 = a^(v,t>), 

for \fv G □ 

Finally, we proceed to derive our GFEM auxiliary space preconditioner. Consider 

aw(u,v) :— / {e{u),e{u))p dx Vu,v£W, 
Jn 

and the mapping U G : W — > V G defined as follows: let u G W be such that 

n 

u = (uk), u k — / ] a k ,j(j>i for k — 1, ■ ■ ■ , d, 
i=i 

then 

2 n 

m g = n G (w) := G V G . 

k=l i=l 

Let denote the space of continuous piecewise quadratic finite elements equipped with 
standard inner product a q (-, •) and define the mapping U q : V G V q by 

U q u G :=u G , for Ug G y gfem . 

We note that n g is the natural inclusion from y GFEM into V q . Associated with a q (-, •), we 
write A q : V q — > (V^)' for its isomorphism. Going back to the original problem ([3]), we define 
the conforming finite element spaces on the mesh Tn- 

V '■= {U'- u = ( u k)l=i, u k G Vi} 



8 
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and 

Vi := jw G V : J u ■ v dx = Wv G Rm| C 

We remark that W can be viewed as VF/RM since Ker(nc) = RM and then W is isomorphic 
to Vj_. To see this, define the operator Qrm : W — > W as follows: for v G W, 

(Qrmv, w) ,u = (v, w) 0) u, w G RM C w. 

Using Qrm, we can regard W as Range(7 — Qrm) where / is an identity operator and Range 
denotes the range of an operator and, thus, W is isomorphic to Vj_. 

The auxiliary space lemma (See [T5l ITS] for a detailed proof) for the linear elasticity 
problem reads: 

Lemma 4.3 (Auxiliary Space Lemma). Assume that H : V q — > W is a surjective and bounded 
linear operator, namely, there exists a positive constant c\ such that Vi>q G V q , 

a w (Uv q ,]Iv q ) < cia q {v q ,v q ). 

Also, we suppose that there exists a positive constant cq such that Vu G W , there is v q G V q 
so that 

v = Uv q and a q (v q ,v q ) < coaiy(t>, t>) 

Then 

(11) Cq 2 aw(u, u) < awCnA' 1 !!* A w u, u) < c\aw{u,u), Vu G W. 

We recall that aw(-, •) < a LE(v,v) <a a>w(', •)■ Also, we see that W and Vj_ are identical. 

For these reasons, it suffices to construct a preconditioner for aw(-, •) on the space W. Let 
n g be an inclusion from V G into V q . We define its Hilbert adjoint operator U q : V q — > V G by 

a q (Il q U G , V q ) = a G (u G , itfqVq), Wu G G V G , V q G V q . 

We note that 

(i G (n G u, n G u) ^ aw(u,u), Vw G w, 

and 

a G (v G , v G ) ^ awCn^va, IL^V?), v G G V G , 

since U G is bijective. Also, the Hilbert adjoint operator (IT^ 1 )''' of II^ 1 can be defined as 
above. Associated with a LE (-,-) in (j3]), we write its isomorphism A LE : V± — > (V±)'. We 
summarize the relations between these spaces and their associated interpolants and dual 
space in the commutative diagram provided in Figure HI 

Let II := n^nj. Then II is surjective operator since W is surjective. Letting v q G V q , we 
have 
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v„ 



■» (v q y 



(nj)* 



V G y (V G Y 

A G 



(n^ 1 )* 



W y W 

A w 



v ± (v ± y 

Ale 



Figure 1. Commutative diagram for Ale and our auxiliary spaces. 



a w (Uv q ,Uvg) = aiy(II G 1 nJi;„n G 1 nJi;g) 

< a G (nJw g , U j v q ) = a q (Jl q U\v q , v q 



< a q (U q Ulv q , Ii q U\v q ) 1/2 a q {v q , v,) 1 ' 2 
= a G (Ulv q ,Ulv q ) 1/2 a q (v q ,v q y/ 2 

< awiU^Ulv^Ulv,) 1 / 2 
= a w (Jlv q ,Ylv q ) 1/2 a q (v q ,v q ) 1/2 . 

Further, let v G W and set v q := U q Il G v G V q . Then, for Wu G W, we have 

a w (fiv q ,u) = a V y(n G 1 nJn (? n G w,M) 
= a G {ulu q u G v, (n- 1 )^) 

= a G (U G v, {IL^yu) 
= a G (v,u). 



Therefore, v = Tlv q . Moreover, 



a q (v q ,v q ) = a q (n q n G v,Il q Il G v) 
= a G (U G v,Il G v) 
< a w (v,v). 

Consequently, we obtain 

(12) ^(n^ntA-^nt)*^ 1 )*^) = ^fi^fiM^) < 1. 



10 
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For notational brevity, let B := HA' 1 !!*. By (JT2J) and Theorem H2J, we get k(BA le ) < A, 
independent of mesh size h. 

We define the norm || • ||i on W as follows: 

ll«lli := + Kl?, 

where u = (111,112) £ W and write for the isomorphism associated with || ■ Wf. Then, using 
the Korn inequality and the Cauchy-Schwarz inequality, we get 

(13) I (e(u),e(u)) F dx^\\u\\l, \/u e W. 

Jn 

These spectral equivalence relations we proposed basically motivate our choice of GFEM- 
based auxiliary space preconditioner. We also remark that the space H 1 ^) has the usual 
norm 

11 11 /11 n2 11 ii2\1/2 

\\u\\H\n) ■= (Nil + \m\L2) ■ 

We describe an algorithm under which several numerical experiments in next subsec- 
tion will be performed. Applying the shape functions on FEM into the bilinear forms 
o>le(', ■)> a w{', ") an d ai(-,-), we obtain the Galerkin matrix A, which represents matrices 
Ale, A w and A\. We consider only a linear system 

Ax = f, 

where / is appropriately chosen. This simplification makes sense by the spectral equiva- 
lent argument we established. The algorithm follows the preconditioned conjugate gradient 
methods in [TB] . 

Algorithm 4.4. 

Compute r = f — Ax , z = 115^11 r 0; and p = z 

For j = 0, 1, until convergence Do: 

r .i+i = <\, - a J A Pj 

z j+1 = IlB q Il r j+ i 

13 j = (r j+u z j+1 )/(r jjZj ) 

Pj+l = Z 3+l + 

EndDo 

Here, II is the matrix representation of II and B q is an approximate inverse of A q where A q is 
the Galerkin matrix arising from the piecewise quadratic FEM. The numerical experiments in 
next section are performed using the stopping criterion jj^J- < 10~ 8 , where is the residual 

of fcth iteration and initial guess Xq = (1, —1, . . . , 1, — 1) T . Taking an account into eigenvalues 
of the preconditioned system, we can estimate the condition number from parameters in the 
conjugate gradient algorithm (see [H] for more detail). 

5. Numerical experiments 

Here we report computed estimates of the condition number of our auxiliary-space pre- 
conditioned system matrix for the linear elasticity model discretized using various choices 
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of mesh spacing on unit square Q = [0, l] 2 . Here, h represents the lengths of the horizontal 
and vertical sides of triangles in the meshes. 




(a) h = i 



(b) h=\ 



{c)h=\ 



(d) h = X 



(e) h = X 



Figure 2. Meshes on the unit square 

Experiment 1. We consider a pair of poisson equations on Q with homogeneous Neumann 
boundaries: 



(14) 



A 
A 



d_ 

dn 







o 

dn , 



u = in Q, 



u = on dQ. 



We note that A\ is the isomorphism associated with the weak formulation of (I14j) on the 
space W . The results in Table [1] implies that the condition number of the preconditioned 
system of A\ is bounded independent of mesh size h. 



h 


i 


i 


i 


i 


i 


4 


s 


if; 


32 




iter 


20 


24 


27 


27 


27 


K 


1.92e+l 


2.80e+l 


3.51e+l 


4.02e+l 


4.39e+l 



Table 1. Condition numbers ^n^nM"^^^ 1 )*^ 



Experiment 2. We deal with a special case when A = 0. Namely, the linear elasticity 
reads 



— div(e(u)) = / in f2 



(15) 



e{u)n 



0, 1< i < 4, 



4. Also, without loss of 

is 



where Tt represent four sides of the unit square for i = 1, . 

generality, we suppose that a — ~, here and later. Assume that the body force / = ( }) i 
defined by 

(16) /i = and f 2 = 0. 

Then the exact solution u G H 2 ± (Q) is 



u 



12 



JAMES BRANNICK AND DURKBIN CHO 



We see that A\y is the isomorphism associated with the weak formulation of ( 115]) on W. In 
Table [21 k — C/c, where C and c are the smallest and largest constants satisfying: 

ca w (u,u) < (BqU.qU.au, n g Il G 'u) < Ca w (u,u),\/u G V±. 
In Table [2j we observe that the preconditioner we devised is an efficient preconditioner for 



h 


1 

4 


i 

s 


i 

ifi 


i 

.32 


i 

64 


iter 


19 


19 


17 


16 


15 


K 


6.47e+0 


6.70e+0 


6.71e+0 


6.67e+0 


6.48e+0 



Table 2. Condition numbers niU^WA- 1 (Ul)* (Uq 1 )* A w ) . 



the special case. Namely, k(BAw) < 1. 

Experiment 3. In this experiment, we study more general linear elasticity problems with 
various A: 

— div(cr(n)) =0 in Cl = unit square 

(17) 

a(u)n =0 1 < I < 4, 

~ ~ r £ 

where cr(u) = e(u) + Xtr(e(u))5 and are defined as in Experiment 2. In Tables [HI |4] and [SI 
k — C/c, where C and c are the smallest and largest constants satisfying: 
ca LE {u,u) < (BqHqUau, U q Uau) < Ca LE {u,u),\/u e V±. 
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24 


1.33e+l 


30 


4.60e+l 


33 


8.60e+l 


-f- 


28 


1.39e+l 


39 


4.51e+l 


50 


8.60e+l 




25 


1.40e+l 


43 


4.57e+l 


55 


8.55e+l 




23 


1.40e+l 


40 


4.56e+l 


53 


8.55e+l 


fi4 


21 


1.39e+l 


36 


4.54e+l 


46 


8.49e+l 



Table 3. Condition numbers K^^njAr 1 ^)*^ 1 )*^^ 



Experiment 4. For small A, we experiment the linear elasticity with homogeneous pure 
traction boundary condition: 

-djy(e(«) + Atr(e(w))5) = / in Q 

(18) 



(e(u) + Xtr(e(u))5)n 



Of 1< £ < 4, 



where Te are also defined as in Experiment 2 and the body force / 



is such that 



u 



x(l — x)y 2 {l — y) 2 sin nx — 



x 2 (l — x) 2 y 2 (l — y) 2 cos Try ' 
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is the exact solution of (I18p . The results reported in Tables H] indicate that the condition 
number of our preconditioned operator remains bounded independent of mesh size and, 
hence, provides for an optimal order preconditioner for our linear elasticity model problem. 
They coincide with our theoretical result that k(BAle) < A, independent of h. We note that 
although our experiments have been carried out on a structured grid, our analysis carries 
over to unstructured simplicial meshes and our scheme is thus expected to perform similar 
for unstructured grids as well. 



h 


A = 1 


A = 5 


A = 10 


iter 


K 


H 1 norm 


iter 


K 


H l norm 


iter 


K 


H L norm 


1 

4 


26 


1.42e+l 


7.71e-3 


32 


4.60c+l 


1.10e-2 


35 


8.60e+l 


1.25e-2 


1 

-+- 


31 


1.39e+l 


4.07e-3 


46 


4.60e+l 


7.56e-3 


56 


8.60e+l 


9.62e-3 


- J f- 


31 


1.41e+l 


1.41e-3 


53 


4.58e+l 


3.48e-3 


68 


8.56e+l 


5.21e-3 


^_ 


31 


1.41e+l 


3.91e-4 


54 


4.59e+l 


1.12e-3 


73 


8.58e+l 


1.92e-3 


64 


31 


1.41e+l 


1.01e-4 


54 


4.59e+l 


3.04e-4 


72 


8.58e+l 


5.49e-4 



Table 4. Condition numbers ^U^U^A^^^^^YAle) and the discrete 
H 1 errors between the exact solution and the numerical solution. 



Challenging experiment. It is well known that the performance of the piecewise linear 
finite elements deteriorates as A approaches oo. In the elasticity literature, it is called the 
locking phenomenon (see [HI H] for more information). To overcome the effects of the locking, 
several methods have been suggested in (2j El QUI El EEl HH [17]. Both our theoretical and 
numerical results show that the GFEM-based auxiliary preconditioner for the linear elasticity 
problem with small A works very efficiently. It would also be worthwhile to investigate 
the numerical results when the preconditioner is applied to linear elasticity for large A. 
We observe in Table |5] that even when A is large enough, the condition number of the 
preconditioned systems is bounded independent of mesh size h and the discretization error 
estimates look reasonable if h is sufficiently small. 



h 


A = 50 (u = 0.49505) 


A = 100 (u = 0.49751) 


A = 500 (u = 0.49950) 


A = 1000 (u = 0.49975) 


iter 


ft 


H L norm 


iter 


ft 


H L norm 


iter 


ft 


H L norm 


iter 


ft 


H L norm 


i 


92 


4.06e+2 


1.39c-2 


109 


8.06e+2 


1.51e-2 


134 


4.01c+3 


1.65e-2 


140 


8.01c+3 


1.66e-2 


-f- 


114 


4.06e+2 


1.04e-2 


144 


8.06e+2 


1.25e-2 


224 


4.01c+3 


1.57e-2 


258 


8.01c+3 


1.64e-2 




135 


4.06e+2 


5.75e-3 


170 


8.06e+2 


8.06e-3 


287 


4.01c+3 


1.31e-2 


365 


8.01c+3 


1.47e-2 




145 


4.05c+2 


2.19c-3 


183 


8.06e+2 


3.69c-3 


322 


4.01e+3 


8.78c-3 


404 


8.01e+3 


1.12e-2 


128 


143 


4.05c+2 


6.42e-4 


192 


8.05e+2 


1.21c-3 


367 


4.00e+3 


4.32c-3 


455 


8.01e+3 


6.36e-3 



Table 5. Condition numbers K(U G 1 U^A~ 1 (W q )*(Jl G 1 )*A LE ) and the discrete 
H 1 errors between the exact solution and the numerical solution. 



6. Concluding remarks 

In this paper, using the auxiliary space technique we have designed a preconditioner for 
the solution of the problem of linear elasticity, which is also based on generalized finite 
element methods. We have proved that for arbitrarily fixed A, the GFEM-based auxiliary 
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preconditioner always works optimally for the system of linear elasticity discretized using 
the lowest finite element shape functions. 

References 

[1] Douglas Arnold. Lecture note on linear plate theory: modeling, analysis, and computation. 

[2] Douglas N. Arnold, Franco Brezzi, Richard S. Falk, and L. Donatella Marini. Locking-free Reissner- 

Mindlin elements without reduced integration. Comput. Methods Appl. Mech. Engrg., 196(37-40) :3660- 

3671, 2007. 

[3] Douglas N. Arnold and Ragnar Winther. Mixed finite elements for elasticity. Numer. Math., 92(3):401- 
419, 2002. 

[4] Ivo Babuska and Manil Suri. Locking effects in the finite element approximation of elasticity problems. 

Numer. Math., 62(4):439-463, 1992. 
[5] Ivo Babuska and Manil Suri. On locking and robustness in the finite clement method. SIAM J. Numer. 

Anal, 29(5):1261-1293, 1992. 
[6] Susanne C. Brenner. A nonconforming mixed multigrid method for the pure traction problem in planar 

linear elasticity. Math. Comp., 63(208) :435-460, S1-S5, 1994. 
[7] Susanne C. Brenner and L. Ridgway Scott. The mathematical theory of finite element methods, volume 15 

of Texts in Applied Mathematics. Springer- Verlag, New York, 1994. 
[8] Susanne C. Brenner and Li-Yeng Sung. Linear finite element methods for planar linear elasticity. Math. 

Comp., 59(200):321-338, 1992. 
[9] Z. Cai, C.-O. Lee, T. A. Manteuffel, and S. F. McCormick. First-order system least squares for lin- 
ear elasticity: numerical results. SIAM J. Sci. Comput., 21(5):1706-1727 (electronic), 2000. Iterative 
methods for solving systems of algebraic equations (Copper Mountain, CO, 1998). 

[10] Zhiqiang Cai, Thomas A. Manteuffel, Stephen F. McCormick, and Seymour V. Parter. First-order 
system least squares (FOSLS) for planar linear elasticity: pure traction problem. SIAM J. Numer. 
Anal, 35(l):320-335 (electronic), 1998. 

[11] D. Cho and L. Zikatanov. A multilevel preconditioning for generalized finite clement method problems 
on unstructured simplicial meshes. J. Numer. Math., 15(3):163-180, 2007. 

[12] Richard S. Falk. Nonconforming finite element methods for the equations of linear elasticity. Math. 
Comp., 57(196):529-550, 1991. 

[13] Leopoldo P. Franca and Rolf Stenberg. Error analysis of Galcrkin least squares methods for the elasticity 
equations. SIAM J. Numer. Anal, 28(6):1680-1697, 1991. 

[14] Chang-Ock Lee. A conforming mixed finite element method for the pure traction problem of linear 
elasticity. Appl. Math. Comput, 93(l):ll-29, 1998. 

[15] S.V. Nepomnyaschikh. Decomposition and fictitious domains methods for elliptic boundary value prob- 
lems. SIAM, Philadelphia, pages 581-613, 1992. Fifth International Symposium on Domain Decompo- 
sition Methods for Partial Differential Equations. 

[16] Yousef Saad. Iterative methods for sparse linear systems. Society for Industrial and Applied Mathemat- 
ics, Philadelphia, PA, second edition, 2003. 

[17] Joachim Schoberl. Multigrid methods for a parameter dependent problem in primal variables. Numer. 
Math., 84(1):97-119, 1999. 

[18] J. Xu. The auxiliary space method and optimal multigrid preconditioning techniques for unstructured 
grids. Computing, 56(3):215-235, 1996. International GAMM- Workshop on Multi-level Methods (Meis- 
dorf, 1994). 

[19] Jinchao Xu. Iterative methods by space decomposition and subspace correction. SIAM Rev., 34(4):581- 
613, 1992. 

Department of Mathematics, The Pennsylvania State University, University Park, PA 
16802, USA. 

E-mail address: brannick@math.psu.edu 

Department of Mathematics, The Pennsylvania State University, University Park, PA 
16802, USA. 

E-mail address: cho@math.psu.edu 



